library("dplyr")
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library("lubridate")
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
library("MASS")
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
library("MuMIn")
library("forecast")
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
source("plot.roots.R")
data <- read.csv("unemployment.csv", header = TRUE)
data <- data[data$State.Area == "California", ]
data <- data %>%
mutate(Total.Unemployment.in.State.Area = as.numeric(gsub("[^0-9]", "", Total.Unemployment.in.State.Area)))
#Subset training to till 2018, testing till 2019
data_training <- data %>%
filter(Year >= 2010 & Year <= 2018)
data_testing <- data %>%
filter(Year >= 2010 & Year <= 2019)
#Make time series
data_training_ts = ts(data_training$Total.Unemployment.in.State.Area, start = c(2010, 1), end = c(2018, 1), frequency = 12)
data_testing_ts = ts(data_testing$Total.Unemployment.in.State.Area, start = c(2010, 1), end = c(2019, 1), frequency = 12)
#Plot time series
plot(data_training_ts, xlab = "Year", ylab = "Unemployment", main = "Figure 1: Unemployment in California from 2010 to 2018", ylim = c(min(data_training_ts), max(data_training_ts)))
#Let's transform data to make it more normal
#Box-Cox Transformation
bcTransform <- boxcox(data_training_ts ~ as.numeric(1:length(data_training_ts)))
bcTransform$x[which(bcTransform$y == max(bcTransform$y))]
## [1] 0.3838384
lambda = bcTransform$x[which(bcTransform$y == max(bcTransform$y))]
data_ts.bc = (1/lambda)*(data_training_ts^lambda-1)
plot(data_ts.bc, main = "Box-Cox Transformation of Unemployment")
hist(data_ts.bc, col = "lightblue", xlab = "Year", ylab = "Unemployment", main = "Box-Cox Transformation Histogram")
#Log Transformation
data_ts.log <- log(data_training_ts)
plot(data_ts.log, main = "Log Transformation of Unemployment")
hist(data_ts.log, col="lightblue", xlab = "Year", ylab = "Unemployment", main = "Log Transformation Histogram")
#Sqrt Transformation
data_ts.sqrt <- sqrt(data_training_ts)
plot(data_ts.sqrt, main = "SQRT Transformation of Unemployment")
hist(data_ts.sqrt, col="lightblue", xlab = "Year", ylab = "Unemployment", main = "Figure 2: SQRT Transformation Histogram")
#The Box-Cox command gives us lambda = 0.3838384. Lambda = 0 corresponds to a Log transformation, while lambda = 0.5 corresponds to a SQRT transformation. Since 0.5 is closer to 0.3838384 and is within the 95% confidence interval, we choose the SQRT transformation. This is also supported by a more symmetric looking histogram.
#Difference to get rid of trend (lag 1)
diff_lag_1_data_ts.transformed = diff(data_ts.sqrt, 1)
plot(diff_lag_1_data_ts.transformed, xlab = "Year", ylab = "Unemployment", main ="De-Trended Unemployment")
abline(h = mean(diff_lag_1_data_ts.transformed),lty = 2)
var(diff_lag_1_data_ts.transformed)
## [1] 13.35966
#Difference to get rid of seasonality (lag 1 + lag 12)
diff_lag_1_and_12_data_ts.transformed = diff(diff_lag_1_data_ts.transformed, 12)
plot(diff_lag_1_and_12_data_ts.transformed, xlab = "Year", ylab = "Unemployment", main ="De-Trended and De-seasonalized Unemployment")
abline(h = mean(diff_lag_1_and_12_data_ts.transformed),lty = 2)
var(diff_lag_1_and_12_data_ts.transformed)
## [1] 23.9294
#Differencing at lag 1 gives us a variance of 3.822456e-05, while differencing at lag 1 and lag 12 gives us a variance of 6.681259e-05. Since the former is a lower variance, we choose to de-trend only.
#De-Trended ACF and PACF
#ACF
acf(diff_lag_1_data_ts.transformed, xlim = c(0, 5), main = " Figure 3: ACF of Unemployment in California from 2010 to 2018", col = "orchid", xlab = "Lag")
#PACF
pacf(diff_lag_1_data_ts.transformed, xlim = c(0, 5), main = "Figure 4: PACF of Unemployment in California from 2010 to 2018", col = "orchid", xlab = "Lag")
#ARIMA looks like the most appropriate model. p is most likely either 2 or 4 since there is a strong AR component, d is either 1 or 2 since we differenced at 1 to remove trend but might need to account for more differencing, and q is definitely 0 since there is no MA component (seasonality).
#ARIMA testing
arima(data_ts.sqrt, order=c(2,1,0), method="ML")
##
## Call:
## arima(x = data_ts.sqrt, order = c(2, 1, 0), method = "ML")
##
## Coefficients:
## ar1 ar2
## 1.6937 -0.7366
## s.e. 0.0687 0.0690
##
## sigma^2 estimated as 1.144: log likelihood = -144.97, aic = 295.93
arima(data_ts.sqrt, order=c(4,1,0), method="ML")
##
## Call:
## arima(x = data_ts.sqrt, order = c(4, 1, 0), method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 ar4
## 2.3097 -2.1781 0.9204 -0.0640
## s.e. 0.1055 0.2472 0.2479 0.1063
##
## sigma^2 estimated as 0.4523: log likelihood = -101.74, aic = 213.48
arima(data_ts.sqrt, order=c(2,2,0), method="ML")
##
## Call:
## arima(x = data_ts.sqrt, order = c(2, 2, 0), method = "ML")
##
## Coefficients:
## ar1 ar2
## 1.2784 -0.7952
## s.e. 0.0619 0.0613
##
## sigma^2 estimated as 0.4604: log likelihood = -99.31, aic = 204.62
arima(data_ts.sqrt, order=c(4,2,0), method="ML")
##
## Call:
## arima(x = data_ts.sqrt, order = c(4, 2, 0), method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 ar4
## 1.3482 -1.1001 0.3960 -0.2519
## s.e. 0.1040 0.1747 0.1759 0.1046
##
## sigma^2 estimated as 0.4315: log likelihood = -96.39, aic = 202.78
AICc(arima(data_ts.sqrt, order=c(2,1,0), method="ML"))
## [1] 296.1931
AICc(arima(data_ts.sqrt, order=c(4,1,0), method="ML"))
## [1] 214.1453
AICc(arima(data_ts.sqrt, order=c(2,2,0), method="ML"))
## [1] 204.8874
AICc(arima(data_ts.sqrt, order=c(4,2,0), method="ML"))
## [1] 203.4565
#The AICc values for the models are as follows: ARIMA(2,1,0) = 296.1931, ARIMA(4,1,0) = 214.1453, ARIMA(2,2,0) = 204.8874, and ARIMA(4,2,0) = 203.4565. Since ARIMA(4,2,0) has the lowest AICc value, we select this model.
#Let's check for stationarity since we know it is invertible (definition of AR)
#ARIMA(2,1,0)
coeff1 <- c(1, -1.6937, 0.7366)
roots1 <- polyroot(coeff1)
plot.roots(NULL, roots1, main = "Roots of AR Polynomial for ARIMA(2,1,0)")
#roots are outside the unit circle, so our model is stationary
#ARIMA(4,1,0)
coeff2 <- c(1, -2.3097, 2.1781, -0.9204, 0.0640)
roots2 <- polyroot(coeff2)
plot.roots(NULL, roots2, main = "Roots of AR Polynomial for ARIMA(4,1,0)")
#roots are outside the unit circle, so our model is stationary
#ARIMA(2,2,0)
coeff3 <- c(1, -1.2784, 0.7952)
roots3 <- polyroot(coeff3)
plot.roots(NULL, roots3, main = "Roots of AR Polynomial for ARIMA(2,2,0)")
#roots are outside the unit circle, so our model is stationary
#ARIMA(4,2,0)
coeff4 <- c(1, -1.3482, 1.1001, -0.3960, 0.2519)
roots4 <- polyroot(coeff4)
plot.roots(NULL, roots4, main = "Roots of AR Polynomial for ARIMA(4,2,0)")
#roots are outside the unit circle, so our model is stationary
#Let's fit the model ARIMA(2,1,0) and get residuals
fit <- arima(data_ts.sqrt, order=c(2,1,0), method="ML")
res1 <- residuals(fit)
#Let's look at the residuals closely
hist(res1, density=20, breaks=20, col="pink", prob=TRUE)
mean <- mean(res1)
standard_deviation <- sqrt(var(res1))
curve(dnorm(x, mean, standard_deviation), add=TRUE )
plot.ts(res1, main = "ARIMA(2,1,0) Residuals")
#Q-Q Plot
fitt <- lm(res1 ~ as.numeric(1:length(res1))); abline(fitt, col="purple")
abline(h=mean(res1), col="lightblue")
qqnorm(res1,main= "Q-Q Plot")
qqline(res1,col="darkblue")
#Seems to follow normal distribution.
#ACF and PACF
acf(res1, lag.max=40)
pacf(res1, lag.max=40)
acf(res1^2, lag.max=40)
#Shapiro-Wilk Test
shapiro.test(res1)
##
## Shapiro-Wilk normality test
##
## data: res1
## W = 0.98112, p-value = 0.1771
#Since p-value = 0.1771, we fail to reject the null hypothesis that the residuals are normally distributed. Therefore the residuals are normally distributed.
#Box-Pierce Test
Box.test(res1, lag = 12, type = c("Box-Pierce"), fitdf = 2)
##
## Box-Pierce test
##
## data: res1
## X-squared = 113.93, df = 10, p-value < 2.2e-16
#Since p-value < 2.2e-16, we reject the null hypothesis that the residuals are white noise. Therefore the residuals are not white noise and have some autocorrelation.
#Ljung-Box Test
Box.test(res1, lag = 12, type = c("Ljung-Box"), fitdf = 2)
##
## Box-Ljung test
##
## data: res1
## X-squared = 122.6, df = 10, p-value < 2.2e-16
#Since p-value < 2.2e-16, we reject the null hypothesis that the residuals are independently distributed. Therefore the residuals are not independent and there is significant autocorrelation.
#McLeod-Li Test
Box.test(res1^2, lag = 12, type = c("Ljung-Box"), fitdf = 0)
##
## Box-Ljung test
##
## data: res1^2
## X-squared = 11.826, df = 12, p-value = 0.4597
#Since p-value = 0.4597, we fail to reject the null hypothesis that the variance of the residuals is white noise. Therefore there is no significant autocorrelation in the variance.
#Yule-Walker Test
ar(res1, aic = TRUE, order.max = NULL, method = c("yule-walker"))
##
## Call:
## ar(x = res1, aic = TRUE, order.max = NULL, method = c("yule-walker"))
##
## Coefficients:
## 1 2 3 4 5 6 7 8
## 0.5198 -0.4311 -0.1273 -0.3891 0.0915 -0.1138 -0.0702 -0.0515
## 9 10 11 12 13
## 0.2757 0.0337 0.0415 -0.0844 0.2884
##
## Order selected 13 sigma^2 estimated as 0.4928
#sigma^2 = 0.4928.
#This model fails the Box-Pierce Test and Ljung-Box Test.
#Let's fit the model ARIMA(4,1,0) and get residuals
fit <- arima(data_ts.sqrt, order=c(4,1,0), method="ML")
res2 <- residuals(fit)
#Let's look at the residuals closely
hist(res2, density=20, breaks=20, col="pink", prob=TRUE)
mean <- mean(res2)
standard_deviation <- sqrt(var(res2))
curve(dnorm(x, mean, standard_deviation), add=TRUE )
plot.ts(res2, main = "ARIMA(4,1,0) Residuals")
#Q-Q Plot
fitt <- lm(res2 ~ as.numeric(1:length(res2))); abline(fitt, col="purple")
abline(h=mean(res2), col="lightblue")
qqnorm(res2,main= "Q-Q Plot")
qqline(res2,col="darkblue")
#Seems to follow normal distribution.
#ACF and PACF
acf(res2, lag.max=40)
pacf(res2, lag.max=40)
acf(res2^2, lag.max=40)
#Shapiro-Wilk Test
shapiro.test(res2)
##
## Shapiro-Wilk normality test
##
## data: res2
## W = 0.9904, p-value = 0.7159
#Since p-value = 0.7159, we fail to reject the null hypothesis that the residuals are normally distributed. Therefore the residuals are normally distributed.
#Box-Pierce Test
Box.test(res2, lag = 12, type = c("Box-Pierce"), fitdf = 2)
##
## Box-Pierce test
##
## data: res2
## X-squared = 25.517, df = 10, p-value = 0.004448
#Since p-value = 0.004448, we reject the null hypothesis that the residuals are white noise. Therefore the residuals are not white noise and have autocorrelation.
#Ljung-Box Test
Box.test(res2, lag = 12, type = c("Ljung-Box"), fitdf = 2)
##
## Box-Ljung test
##
## data: res2
## X-squared = 28.573, df = 10, p-value = 0.00146
#Since p-value = 0.00146, we reject the null hypothesis that the residuals are independently distributed. Therefore the residuals are not independent and there is significant autocorrelation.
#McLeod-Li Test
Box.test(res2^2, lag = 12, type = c("Ljung-Box"), fitdf = 0)
##
## Box-Ljung test
##
## data: res2^2
## X-squared = 19.714, df = 12, p-value = 0.0727
#Since p-value = 0.0727, we fail to reject the null hypothesis that the variance of the residuals is white noise. Therefore there is no significant autocorrelation in the variance.
#Yule-Walker Test
ar(res2, aic = TRUE, order.max = NULL, method = c("yule-walker"))
##
## Call:
## ar(x = res2, aic = TRUE, order.max = NULL, method = c("yule-walker"))
##
##
## Order selected 0 sigma^2 estimated as 0.4669
#sigma^2 = 0.4669.
#This model fails the Box-Pierce Test and Ljung-Box Test.
#Let's fit the model ARIMA (2,2,0) and get residuals
fit <- arima(data_ts.sqrt, order=c(2,2,0), method="ML")
res3 <- residuals(fit)
#Let's look at the residuals closely
hist(res3, density=20, breaks=20, col="pink", prob=TRUE)
mean <- mean(res3)
standard_deviation <- sqrt(var(res3))
curve(dnorm(x, mean, standard_deviation), add=TRUE )
plot.ts(res3, main = "ARIMA(2,2,0) Residuals")
#Q-Q Plot
fitt <- lm(res3 ~ as.numeric(1:length(res3))); abline(fitt, col="purple")
abline(h=mean(res3), col="lightblue")
qqnorm(res3,main= "Q-Q Plot")
qqline(res3,col="darkblue")
#Seems to follow normal distribution.
#ACF and PACF
acf(res3, lag.max=40)
pacf(res3, lag.max=40)
acf(res3^2, lag.max=40)
#Shapiro-Wilk Test
shapiro.test(res3)
##
## Shapiro-Wilk normality test
##
## data: res3
## W = 0.99178, p-value = 0.82
#Since p-value = 0.82, we fail to reject the null hypothesis that the residuals are normally distributed. Therefore the residuals are normally distributed.
#Box-Pierce Test
Box.test(res3, lag = 12, type = c("Box-Pierce"), fitdf = 2)
##
## Box-Pierce test
##
## data: res3
## X-squared = 22.382, df = 10, p-value = 0.01327
#Since p-value = 0.01327 we reject the null hypothesis that the residuals are white noise. Therefore the residuals are not white noise and have some autocorrelation.
#Ljung-Box Test
Box.test(res3, lag = 12, type = c("Ljung-Box"), fitdf = 2)
##
## Box-Ljung test
##
## data: res3
## X-squared = 25.184, df = 10, p-value = 0.005008
#Since p-value = 0.005008, we reject the null hypothesis that the residuals are independently distributed. Therefore the residuals are not independent and there is significant autocorrelation.
#McLeod-Li Test
Box.test(res3^2, lag = 12, type = c("Ljung-Box"), fitdf = 0)
##
## Box-Ljung test
##
## data: res3^2
## X-squared = 13.144, df = 12, p-value = 0.3587
#Since p-value = 0.3587, we fail to reject the null hypothesis that the variance of the residuals is white noise. Therefore there is no significant autocorrelation in the variance.
#Yule-Walker Test
ar(res3, aic = TRUE, order.max = NULL, method = c("yule-walker"))
##
## Call:
## ar(x = res3, aic = TRUE, order.max = NULL, method = c("yule-walker"))
##
##
## Order selected 0 sigma^2 estimated as 0.4993
#sigma^2 = 0.4993
#This model fails the Box-Pierce Test and Ljung-Box Test.
#Let's fit the model ARIMA(4,2,0) and get residuals
fit <- arima(data_ts.sqrt, order=c(4,2,0), method="ML")
res4 <- residuals(fit)
#Let's look at the residuals closely
hist(res4, density=20, breaks=20, col="pink", prob=TRUE)
mean <- mean(res4)
standard_deviation <- sqrt(var(res4))
curve(dnorm(x, mean, standard_deviation), add=TRUE )
plot.ts(res4, main = "ARIMA(4,2,0) Residuals")
#Q-Q Plot
fitt <- lm(res4 ~ as.numeric(1:length(res4))); abline(fitt, col="purple")
abline(h=mean(res4), col="lightblue")
qqnorm(res4,main= "Q-Q Plot")
qqline(res4,col="darkblue")
#Seems to follow normal distribution.
#ACF and PACF
acf(res4, lag.max=40)
pacf(res4, lag.max=40)
acf(res4^2, lag.max=40)
#Shapiro-Wilk Test
shapiro.test(res4)
##
## Shapiro-Wilk normality test
##
## data: res4
## W = 0.98406, p-value = 0.2903
#Since p-value = 0.2903, we fail to reject the null hypothesis that the residuals are normally distributed. Therefore the residuals are normally distributed.
#Box-Pierce Test
Box.test(res4, lag = 12, type = c("Box-Pierce"), fitdf = 2)
##
## Box-Pierce test
##
## data: res4
## X-squared = 11.622, df = 10, p-value = 0.3111
#Since p-value = 0.3111, we fail to reject the null hypothesis that the residuals are white noise. Therefore the residuals are white noise and have no autocorrelation.
#Ljung-Box Test
Box.test(res4, lag = 12, type = c("Ljung-Box"), fitdf = 2)
##
## Box-Ljung test
##
## data: res4
## X-squared = 13.133, df = 10, p-value = 0.2163
#Since p-value = 0.2163, we fail to reject the null hypothesis that the residuals are independently distributed. Therefore the residuals are independent and there is no significant autocorrelation.
#McLeod-Li Test
Box.test(res4^2, lag = 12, type = c("Ljung-Box"), fitdf = 0)
##
## Box-Ljung test
##
## data: res4^2
## X-squared = 10.735, df = 12, p-value = 0.5518
#Since p-value = 0.5518, we fail to reject the null hypothesis that the variance of the residuals is white noise. Therefore there is no significant autocorrelation in the variance.
#Yule-Walker Test
ar(res4, aic = TRUE, order.max = NULL, method = c("yule-walker"))
##
## Call:
## ar(x = res4, aic = TRUE, order.max = NULL, method = c("yule-walker"))
##
##
## Order selected 0 sigma^2 estimated as 0.4698
#sigma^2 = 0.4698
#This model passes all tests.
#Since ARIMA(4,2,0) is the only model that checks all the diagnostic criteria and has the lowest AICc value, we are confident this is the appropriate model.
#Equation for selected model ARIMA(4,2,0)
model <- arima(data_ts.sqrt, order=c(4,2,0), method="ML")
summary(model)
##
## Call:
## arima(x = data_ts.sqrt, order = c(4, 2, 0), method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 ar4
## 1.3482 -1.1001 0.3960 -0.2519
## s.e. 0.1040 0.1747 0.1759 0.1046
##
## sigma^2 estimated as 0.4315: log likelihood = -96.39, aic = 202.78
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.06814014 0.6852723 0.5639344 -0.003994467 0.04696552 0.08781038
## ACF1
## Training set -0.01794226
#Δ_2 sqrt(U_t) = (1 - 1.3482B + 1.1001B^2 - 0.3960B^3 + 0.2519B^4)Z_t
#Sqrt transformed data plot with 12 forecasts
fit <- arima(sqrt(data_training_ts), order=c(4,2,0), method="ML")
forecast(fit)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Feb 2018 913.6224 912.7806 914.4642 912.3349 914.9098
## Mar 2018 909.6700 906.7285 912.6115 905.1713 914.1686
## Apr 2018 905.8627 899.7143 912.0112 896.4596 915.2659
## May 2018 901.3192 891.3253 911.3132 886.0348 916.6036
## Jun 2018 895.6559 881.6620 909.6498 874.2540 917.0577
## Jul 2018 889.0394 871.2290 906.8499 861.8008 916.2781
## Aug 2018 882.0419 860.7158 903.3680 849.4264 914.6573
## Sep 2018 875.3210 850.6907 899.9513 837.6522 912.9898
## Oct 2018 869.2970 841.3616 897.2324 826.5735 912.0205
## Nov 2018 863.9973 832.5220 895.4726 815.8600 912.1346
## Dec 2018 859.1129 823.6968 894.5290 804.9486 913.2771
## Jan 2019 854.1979 814.4007 893.9952 793.3333 915.0626
## Feb 2019 848.8962 804.3659 893.4264 780.7931 916.9993
## Mar 2019 843.0886 793.6285 892.5488 767.4458 918.7314
## Apr 2019 836.9080 782.4587 891.3574 753.6349 920.1811
## May 2019 830.6354 771.2022 890.0685 739.7403 921.5304
## Jun 2019 824.5462 760.1182 888.9742 726.0121 923.0803
## Jul 2019 818.7853 749.2819 888.2886 712.4891 925.0815
## Aug 2019 813.3226 738.5839 888.0613 699.0196 927.6255
## Sep 2019 807.9966 727.8103 888.1830 685.3622 930.6310
## Oct 2019 802.6108 716.7582 888.4633 671.3106 933.9109
## Nov 2019 797.0291 705.3270 888.7312 656.7828 937.2754
## Dec 2019 791.2284 693.5475 888.9094 641.8383 940.6186
## Jan 2020 785.2897 681.5465 889.0330 626.6281 943.9513
pred.tr <- predict(fit, n.ahead = 12)
U.tr= pred.tr$pred + 2*pred.tr$se
L.tr= pred.tr$pred - 2*pred.tr$se
#Full plot
ts.plot(sqrt(data_testing_ts), xlim = c(2010, 2019), ylim=c(750, 1500), main = "Transformed Unemployment with Forecasts")
lines(U.tr, col="lightblue", lty="dashed")
lines(L.tr, col="lightblue", lty="dashed")
points(seq(from = 2018, to = 2019, length.out=12), pred.tr$pred, col="orchid")
#Zoomed plot
ts.plot(sqrt(data_testing_ts), xlim = c(2017, 2019), ylim=c(750, 1500), main = "Figure 5: Transformed Unemployment with Forecasts (Zoomed)")
lines(U.tr, col="lightblue", lty="dashed")
lines(L.tr, col="lightblue", lty="dashed")
points(seq(from = 2018, to = 2019, length.out=12), pred.tr$pred, col="orchid")
#We can clearly see that the transformed data from 2019 falls within the prediction interval even if the forecasts are not exactly on the data.
#Original data plot with 12 forecasts
#Full plot
pred.orig <- pred.tr$pred ^ 2
U = U.tr ^ 2
L = L.tr ^ 2
ts.plot(data_testing_ts, xlim = c(2010, 2019), ylim=c(100000, 2500000), main = "Actual Unemployment with Forecasts")
lines(U, col="lightblue", lty="dashed")
lines(L, col="lightblue", lty="dashed")
points(seq(from = 2018, to = 2019, length.out=12), pred.orig, col="orchid")
#Zoomed plot
pred.orig <- pred.tr$pred ^ 2
U = U.tr ^ 2
L = L.tr ^ 2
ts.plot(data_testing_ts, xlim = c(2017, 2019), ylim=c(600000, 1200000), main = "Figure 6: Actual Unemployment with Forecasts (Zoomed)")
lines(U, col="lightblue", lty="dashed")
lines(L, col="lightblue", lty="dashed")
points(seq(from = 2018, to = 2019, length.out=12), pred.orig, col="orchid")
#We can clearly see that the actual data from 2019 falls within the prediction interval even if the forecasts are not exactly on the data.